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Traveling  wave  modes  of  a  plane  layered  anelastic  earth 


Robert  I.  Odom 

Applied  Physics  Laboratory  and  Department  of  Earth  and  Space  Sciences 
University  of  Washington,  1013  NE  40th  St.,  Seattle,  WA  98105 


Summary.  Incorporation  of  attenuation  into  the  normal  mode  sum  representations  of  seismic 
signals  is  commonly  effected  by  applying  perturbation  theory.  This  is  fine  for  weak  attenuation,  but 
problematic  for  stronger  attenuation.  In  this  work  modes  of  the  anelastic  medium  are  represented 
as  complex  superpositions  of  elastic  eigenfunctions.  For  the  P-SV  system  a  generalized  eigenvalue 
equation  for  the  complex  eigenwavenumbers  and  complex  coefficients  used  to  construct  the  anelastic 
eigenfunctions  is  derived.  The  generalized  eigenvalue  problem  for  the  P-SV  problem  is  exactly 
linear  in  the  eigenwavenumber  at  the  expense  of  doubling  the  dimension.  The  SH  problem  is  exactly 
linear  in  the  square  of  the  eigenwavenumber.  This  is  in  contrast  to  a  similar  standing  wave  problem 
for  the  earth  free  oscillations  (Tromp  andDahlen,  1990).  Attenuation  is  commonly  incorporated  into 
synthetic  seismogram  calculations  by  introduction  of  complex  frequency  dependent  elastic  moduli. 
The  moduli  depend  nonlinearly  on  the  frequency.  The  independent  variable  in  the  standing  wave 
free  oscillation  problem  is  the  frequency,  which  makes  the  eigenvalue  problem  nonlinear.  The 
choice  of  the  wavenumber  as  the  independent  variable  for  the  traveling  wave  problem  leads  to  a 
linear  problem.  The  Earth  model  may  be  transversely  isotropic.  Compressional  waves,  and  both 
polarizations  of  shear  waves  (SV,  SH)  are  treated. 

Key  words:  Seismic  attenuation;  Seismic  anisotropy;  Surface  waves  and  free  oscillations;  Theo¬ 
retical  seismology,  Wave  propagation 
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1.  Introduction 

A  surface- wave  seismic  or  shallow-water  ocean  seismo-acoustic  signal  can  be  conveniently  repre¬ 
sented  as  a  linear  superposition  of  orthogonal  modes  for  moderately  low  frequencies.  Our  definition 
of  “moderately  low”  is  model  dependent,  by  which  we  mean  that  the  model  is  a  “few”  wavelengths 
thick.  A  “few”  is  somewhat  arbitrary,  but  may  be  as  many  as  100.  The  main  point  is  that  the 
frequency  content  of  the  signal  is  not  high  enough  to  justify  using  ray  theory. 

When  modal  representations  of  the  seismic  or  seismo-acoustic  wave-field  have  been  applied 
to  surface-waves  of  the  solid  Earth  or  shallow  water  sound  propagation  by  the  ocean  acoustics 
community,  it  has  been  fairly  common  practice  to  include  the  effects  of  attenuation  as  a  first  order 
perturbation  to  modal  eigenvalues  (e.g.  Ingenito,  1973;  Aki  and  Richards,  1980;  Zhou,  1985).  First 
order  perturbation  theory  ignores  anelastic  coupling  between  modes  and  requires  that  kj  ( Qdk )  <C  1 
where  k  is  the  unperturbed  wavenumber,  8k  is  the  unperturbed  wavenumber  spacing  and  Q  is  the 
spatial  quality  factor.  Because  at  low  frequencies  a  significant  fraction  of  a  shallow  water  seismo- 
acoustic  or  regional  seismic  signal  may  propagate  in  low  Q  bottom  sediments,  or  low  shear  Q  upper 
mantle  k/(Q8k)  can  be  (9(1).  This  makes  the  use  of  perturbation  theory  invalid,  and  can  introduce 
serious  errors  in  mode  sum  acoustic  signal  synthesis.  Ewing  et  al.  (1992)  report  shear  Q  values 
in  the  range  20  -  50  for  continental  shelf  sediments  off  the  New  Jersey.  Lebedev  and  Nolet  (2003) 
found  shear  Q  as  low  as  40  in  the  upper  mantle. 

The  severity  of  the  error  resulting  from  the  improper  treatment  of  attenuation  in  mode  sum 
signal  synthesis  has  been  graphically  illustrated  by  Day  et  al.  (1989).  Day  et  al.  (1989)  calculated 
synthetic  seismograms  for  stratified  models  consisting  of  a  high  Q  layer  over  a  layered  half  space. 
The  shear  Q  of  the  underlying  half  space  was  lower  than  the  Q  in  the  overlying  layer,  and  half  space 
shear  speeds  were  lower  than  the  compressional  wave  speed  in  the  overlying  layer.  The  signals 
calculated  from  modal  summation  dramatically  overestimate  the  effect  of  the  low  shear  Q  on  the 
complete  signal.  The  mode  summation  results  were  compared  with  the  results  from  a  wavenumber 
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integration  routine  (Apsel,  R.J.  and  J.E.  Luco,  1983),  that  is  similar  to  SAFARI/OASES  (Schmidt 
and  Tango,  1986)  in  common  use  in  ocean  acoustics.  Serious  errors  in  the  mode  sum  seismograms 
are  traced  by  Day  et  al.  (1989)  to  the  way  in  which  perturbation  theory  is  applied  to  treat  the  anelastic 
problem.  Specifically,  the  difficulty  occurs  with  the  use  of  the  unperturbed  elastic  eigenfunctions 
in  the  anelastic  problem.  Because  the  problem  is  with  the  eigenfunctions  themselves,  it  cannot 
be  repaired  with  higher  order  perturbation  corrections  to  the  eigenvalues.  The  problem  does  not 
appear  when  using  wavenumber  integration  routines,  because  the  attenuation  is  incorporated  directly 
through  the  use  of  complex,  frequency  dependent  compressional  and  shear  speeds. 

If  one  wishes  to  retain  the  physical  insight  inherent  in  a  modal  representation  of  the  wave-field, 
and  properly  incorporate  the  effects  of  anelasticity,  one  approach  is  to  invoke  the  correspondence 
principle  (e.g.  Leitman  and  Fisher,  1973),  and  solve  for  the  anelastic  modes  directly.  The  corre¬ 
spondence  principle  states  that  the  equations  of  motion  for  a  linear  viscoelastic  material  are  just 
the  equations  for  a  perfectly  elastic  material  with  the  elastic  moduli  replaced  with  the  complex, 
frequency  dependent  anelastic  moduli.  The  anelastic  moduli  must  be  frequency  dependent  and 
satisfy  the  Kramers-Kronig  relations  to  preserve  causality. 

The  correspondence  principle  approach  has  been  used  by  Yuen  and  Peltier  (1982)  and  Buland  et 
al.  (1985)  to  model  aspects  of  the  free  oscillations  of  the  whole  earth.  To  a  limited  extent,  it  has  also 
been  applied  to  shallow  water  propagation  problems.  Bucker  and  Morris  (1965)  employed  the  cor¬ 
respondence  principle  to  solve  for  the  anelastic  eigenwavenumbers,  and  model  the  propagation  loss 
for  a  shallow  water  problem  with  a  fairly  simple  structure.  It  is  standard  procedure  in  wavenumber 
integration  approaches  to  wave-field  modeling,  e.g.  Apsel,  R.J.  and  J.E.  Luco,  (1983),  Schmidt 
and  Tango  (1986)  where  arbitrary  attenuation  is  incorporated  by  introducing  complex,  frequency 
dependent  elastic  moduli. 

There  has  been  work  on  directly  solving  for  the  complex  modes  of  a  plane  layered  fluid-elastic 
medium.  Ivansson  and  Karasalo  (1992,  1993)  and  Ivansson  (1997)  have  published  a  numerical 
algorithm  based  on  the  winding  number  theorem  from  complex  analysis  to  directly  search  for  the 
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poles  of  the  anelastic  modes  in  the  complex  plane.  The  complex  plane  is  tiled  with  boxes,  and 
contour  integrals  are  performed  numerically  around  each  box.  They  determine  the  number  of  poles 
within  each  box,  and  iteratively  refine  the  search  until  each  box  contains  a  single  pole.  A  final 
refinement  is  effected  by  switching  to  polar  coordinates,  and  integrating  around  circular  contours 
to  more  precisely  isolate  the  anelastic  mode  eigenvalues. 

We  adopt  a  different  approach  in  that  we  represent  the  anelastic  modes  as  a  complex  superpo¬ 
sition  of  elastic  eigenfunctions.  The  effects  due  to  anelastic  mode  coupling  are  explicitly  included 
and  there  is  no  restriction  on  the  magnitude  of  the  damping.  Our  approach  is  a  traveling  wave 
adaptation  of  Tromp  and  Dahlen  (1990),  who  derived  an  elegant  solution  for  the  free  oscillations 
of  an  anelastic  spherical  earth  in  terms  of  the  elastic  eigenfunctions  and  eigenfrequencies. 

We  have  derived  a  generalized  eigenvalue  equation  for  the  complex  eigenwavenumbers  and 
complex  coefficients  used  in  the  superposition  of  the  elastic  eigenfunctions  to  construct  the  anelas¬ 
tic  eigenfunctions.  Our  generalized  eigenvalue  equation  is  strictly  linear  for  the  complex  anelastic 
eigenwavenumbers.  This  is  in  contrast  to  the  nonlinear  eigenvalue  equation  for  the  anelastic  eigen¬ 
frequencies  of  the  free  oscillations  of  the  earth  (Dahlen,  1981;  Tromp  and  Dahlen,  1990).  The 
reason  for  this  difference  is  our  choice  of  the  frequency  ft)  as  the  independent  variable  in  the  dis¬ 
persion  relation.  Because  of  the  standing  wave  nature  of  the  earth  free  oscillation  problem,  ft)  is 
taken  as  the  dependent  variable  in  the  dispersion  relation.  Since  the  anelastic  moduli  are  frequency 
dependent,  the  eigenvalue  equation  for  the  anelastic  free  oscillations  is  nonlinear.  Our  derivation 
also  includes  the  effects  of  vertical  transverse  isotropy,  which  has  a  single  vertical  axis  of  symmetry. 
A  particular  feature  of  transversely  isotropic  media  is  that  the  P-SV  motion  still  decouples  from 
the  SH  motion.  This  is  not  true  for  more  general  anisotropy.  The  case  of  transverse  isotropy  is 
particularly  relevant  for  bottom-interacting  shallow  water  sound  propagation.  Berge  et  al.(1991) 
felt  that  additional  azimuthal  anisotropy  induced  by  ripples  in  the  sediment  surface  or  cracks  would 
be  very  weak. 

Bottom-interacting  shallow  water  sound  propagation  requires  special  handling  of  the  fluid- 
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solid  interface  in  the  P-SV  problem.  There  are  no  inherent  difficulties  in  treating  this  interface.  The 
procedure  is  straightforward  and  discussed  in  the  next  section  of  the  manuscript. 

We  include  the  SH  problem  for  completeness,  but  discuss  it  only  briefly.  The  derivation  of  the 
generalized  eigenvalue  equation  for  the  anelastic  eigenwavenumbers  will  be  sketched. 


2.  Definitions  and  notation 


We  adopt  the  notational  conventions  of  Takeuchi  and  Saito  (1972),  who  treat  seismic  surface-waves 
and  free  oscillations  explicitly  for  a  transversely  isotropic  Earth.  Our  coordinate  system  is  a  right 
handed  coordinate  system  with  the  propagation  direction  along  the  x  axis,  y  positive  into  the  paper, 
and  z  positive  upward.  As  mentioned  the  P-SV  (Rayleigh)  motion  decouples  from  the  SH  (Love) 
motion.  The  perfectly  elastic  displacements  «,  and  stresses  Ojj  for  P-SV  are 


UX  = 

Uy  =  0  (1) 

uz  =  yi(z;G>, *)*'■(“»-**> 


and 


&XX 

i  ( at  —  kx ) 

°yy 

■ 

2N)y3^  ei(cot 

®zz 

=  y2ei^m~kx) 

&ZX 

=  -iyAei{m-kx) 

Gyz 

—  &xy  =  0- 

(2) 
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The  boundary  conditions  on  the  y,  for  P-SV  are 
( i  =  1,2,  3,4)  continuous 

yi  =  =  0  at  the  free  surface  z  —  0  (3) 

>’/  (i  =  1,2, 3,4)  -A  0  as  z  4  — 

The  displacements  and  stress  for  SH  motion  are 


Ux  —  u i  —  0 

uy  =  yi  {r,(o,k)e^m~kx)  (4) 


&y  z 

II 

:( at  —  kx) 

&xy 

=  —  ik'Hyxe1^' 

eot  —  kx) 

(5) 

At 

=  Oyy  =  Czz 

=  &ZX  = 

In  addition,  we  introduce  the  definition  for  y2  so  that 


crv-  =  y2ei{a,-kx)  =  L— 

dz 


The  boundary  conditions  for  the  y,  for  SH  are 


(6) 


yi,y2  continuous 

y2  =  0  at  the  free  surface  z  =  0  (7) 

yi,y2  ->■  0  as  z  4  — 


(U  is  the  real  angular  frequency,  p  is  the  density,  k  is  the  horizontal  wavenumber  for  a  perfectly 
elastic  medium,  and  A,  C,  F,  L  and  N  are  the  five  elastic  moduli  in  Love’s  (1944)  notation  necessary 
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to  characterize  a  transversely  isotropic  medium.  When 

A  =  C  =  A+2/i  L  =  N  =  ju  F  =  A  (8) 

the  medium  is  isotropic. 

Fluid  layers  and  the  boundary  between  fluid  and  solid  layers  require  special  treatment  since  fluids 
do  not  support  shear.  In  the  fluid 

A  =  C  =  F  =  A  L  =  N  =  0  (9) 

The  displacements  and  stresses  in  the  fluid  are 

ux  =  -iy3(r,(0,k)ei{cot~kx)  = 

co-p' 

Uy  =  0  (10) 

uz  =  y\{z',(0,k)el^m~kx^ 

and 


6 XX 

=  A  -  Ay3^)  e^w-kx^ 

®zz 

=  yiei((0t-kx) 

(ID 

&ZX 

=  -iy4ei(0}t-kx^ 

ayy 

—  &yz  ~  ®xy  =  0- 

The  boundary  conditions  on  the  y,  for  fluids  are 

yi  (i  —  1,2)  continuous  at  the  fluid  —  solid  boundary 
k 

y3  — - y?  in  the  solid  at  the  fluid  —  solid  boundary  (12) 

C0ZP' 

y4  =0  in  the  fluid  and  at  the  fluid  —  solid  boundary  in  the  solid  . 
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In  order  to  minimize  the  notational  overhead,  we  make  no  distinction  between  the  y,-’s  for  the 
P-SV  and  SH  problems.  We  will  concentrate  mainly  on  the  P-SV  problem.  The  meaning  of  the 
yi  s  will  be  clear  from  context.  The  P-SV  and  SH  problems  are  treated  separately  since  they  do 
not  couple  in  transversely  isotropic  media  with  a  vertical  symmetry  axis.  The  above  equations  of 
motion  for  a  perfectly  elastic  medium  may  be  represented  in  first  order  form  as 

d-btf;L  —  (13) 

where  the  subscripts  R  (Rayleigh),  L  (Love)  indicate  whether  we  are  referring  to  the  P-SV,  Eq.  (1) 
-  (3),  or  the  SH,  Eq.  (4)  -  (7),  systems  of  equations.  The  vectors  b^L  and  the  matrices  Mrl  are 
defined  for  P-SV  motion  as 

bR  =  (yi,j2,y3,j4)T  (14) 

and 


Mr 


o  cr1  k¥/c  o 

~(02p  0  Ok 

-k  0  0  L"1 

0  —kF/C  [k2  (A  —  F2/C)  —  (02p]  0 


for  SH  motion  as 


bL  =  {y\,yi)T 


(15) 


(16) 


and 
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0  L  1 


[k2 N  -  arp ) 


(17) 


and  for  fluids  as 


b  Fluid  =  {y\,yi)T 


(18) 


and 


Mpiuid 


~(o2p  0 


(19) 


From  this  point  on,  we  drop  the  subscripts  R,  L,  and  Fluid  on  the  vector  b  and  the  matrix  M.  It  will 
be  apparent  from  context  which  system  we  mean.  The  following  development  will  be  for  the  P-SV 
system.  Analogous  results  for  the  SH  system  are  summarized  at  the  end  of  the  paper.  The  details 
of  incorporating  fluid-solid  boundaries  are  well  known.  Excellent  treatments  are  in  Takeuchi  and 
Saito  (1972)  and  Aki  and  Richards  (1980,  pp280-281). 


There  are  inherent  symmetries  in  the  equations  of  motion  (Kennett  et  al.,  1978  and  Thomson  et 
al.,  1986)  that  can  be  exploited  to  construct  compact  expressions  useful  for  very  efficiently  deriving 
the  elastic  wave  dispersion  relation,  orthogonality  relationships  and  other  quantities.  Define  the 
matrices  R,  S,  A  and  S  as 


A  0 

i 

o 

M 

and  S  = 

0  A 

0  S 

(20) 
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and 


0  1 

0  1 

A  = 

-1  0 

and 

'h'  — 

— 

1  0 

(21) 


Using  the  matrices  defined  above,  we  can  form  various  compositions  of  two  stress-displacement 
fields.  For  example  for  P-SV,  we  can  form 


<9z(brSb)  -  br  [MrS  +  SM]  b 


(22) 


whereby  we  operate  from  the  left  and  the  right  with  the  same  stress  displacement  vector  b.  The 
perfectly  elastic  modes  are  then  b,„  for  m  =  0, 1,2, The  completeness  of  the  bm  for  the  case 
of  a  plane  layered  medium  with  a  free  surface  and  a  rigid  lower  boundary  has  been  proved  by 
Kirrmann(1995).  The  only  approximation  is  that  in  any  practical  implementation,  we  employ  a 
finite  number  of  modes.  When  both  sides  are  integrated  with  respect  to  z  from  —  °o  to  0.  we  obtain 
the  dispersion  relation  for  Rayleigh  waves  in  a  transversely  isotropic  medium  (Aki  and  Richards, 
1980). 


CD  Ij  —  k  I2  T  kl  3  T- 14 


(23) 


where 


It 


p(y\2  +y32)dz 


(24) 


i2 


(Lyi2  +  A y32)dz 


I3 


I4  = 


T  dy-i 
+  L— ^ 
dz 


dz. 


(25) 

(26) 

(27) 
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The  left  hand  side  of  Eq.(22)  is  a  perfect  differential  in  z,  and  after  integration  it  vanishes  when  the 
boundary  conditions  are  applied. 

3.  Derivation  of  the  generalized  eigenvalue  equation 

The  complex  generalized  eigenvalue  equation  for  the  complex  eigenwavenumbers  is  derived  in 
a  straightforward  manner.  Invoking  the  correspondence  principle,  we  represent  the  equations  of 
motion  for  an  anelastic  transversely  isotropic  medium  as 

dzc  —  Me  (28) 

where  c  and  M  represent  the  stress-displacement  vector  and  wave  operator  matrix,  respectively  for 
anelastic  media.  For  our  definition  of  M,  we  take 

C-1  fcF/C  0 

0  Ok 

(29) 

0  0  F"1 

kF/C  [k2  (A  -  F2/C)  -  a2p]  0 

The  symbols  A,  C,  F,  L,  and  N  are  the  five  complex  frequency  dependent  elastic  moduli  for  an 
anelastic  transversely  isotropic  solid;  K  is  the  eigenwavenumber  for  the  anelastic  solid;  o  is  the 
frequency,  which  we  take  to  be  real  for  propagating  waves. 

We  take 

c  =  cn  (30) 

where  cn  is  an  eigenfunction  of  the  anelastic  medium  and  is  a  solution  of  Eq.  (28).  In  addition,  we 
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represent  c„  as  a  complex  linear  superposition  of  the  eigenfunctions  of  the  perfectly  elastic  problem 

oo 

=  Qmtt^m-  (31) 

m= 0 

Qmn  is  the  matrix  of  complex  coefficients  that  transforms  the  elastic  eigenfunctions  b,„  to  the 
anelastic  solutions. 

Employing  the  matrix  R  defined  above,  we  form  the  composition  of  an  anelastic  mode  c„  as 
represented  by  Eq.  (31)  and  an  elastic  mode  b„ 

dz(blRcn)  =  bft  [MrR  +  RM]  cn,  (32) 

and  integrate  over  z  from  —  oo  to  0.  The  elastic  and  anelastic  problems  satisfy  the  same  boundary 
conditions,  so  the  left  hand  side  of  Eq.  (32)  vanishes  after  the  integration.  By  assuming  that  the 
elastic  eigenfunction  b„  and  the  anelastic  eigenfunction  cn  have  the  same  real  frequency  so  that 

ft)2  =  a2,  (33) 

we  arrive  at  the  following  infinite  generalized  quadratic  eigenvalue  equation  for  the  anelastic  eigen- 
wavenumbers  Kn 


Aq„  +  Kn  B  q„  +  k2  Cq„  =  0 


(34) 


where 


A  =  -  f°  (A  —  F2/C)  y3(w)y3(w)| 

J  — OO  V  s 
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+  kn  (yi("Wm)  +  y4(n)yiw)  -I(y2{n)y3{m)  +  y3{n)y2{m))  , 

-  (C"1  -C-1)y2{n)y2{m)  ~  (L_1  -  L-1)y4(BWm)}  (35) 

b  =  J°  {y^nWm)+y4{n)yi{m))  -  E(j2(n)y3(w)+y3(M)y2(w))  (36) 

C  =  /°  [(A  -  F2/C) b(m)B(w)1  dz.  (37) 

J  — oo  L 

As  Eq.  (38)  shows,  the  eigenvectors  q„  are  the  columns  of  the  matrix  Q 

Q  -  (-..,  qM,  •••)•  (38) 

By  making  the  assignment 

IrM  —  KnIqn,  (39) 


the  quadratic  generalized  eigenvalue  problem  Eq.  (34)  can  be  converted  to  a  linear  generalized 
eigenvalue  problem  (Garbow  et  al.,  pp. 49-50,  1977)  at  the  expense  of  doubling  the  dimension  of 
the  system 

(40) 

Equation  (40)  is  the  main  result  of  this  note.  The  solution  of  this  linear  generalized  matrix  eigen¬ 
value  problem  yields  the  complex  eigenwavenumbers  Kn  for  the  modes  of  an  anelastic  transversely 
isotropic  medium  and  the  eigenvectors  qn.  As  Eq.(38)  shows,  the  eigenvectors  qn  of  Eq.  (40)  are 
the  columns  of  the  transformation  matrix  Qmn  used  to  construct  the  anelastic  eigenfunctions  from 
the  elastic  eigenfunctions  from  Eq.  (31).  The  linearity  of  Eq.  (40)  is  an  important  point  and  should 
be  contrasted  with  the  result  derived  by  Tromp  and  Dahlen  (1990)  for  the  free  oscillations  of  the 
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earth.  Their  equation  (3.5)  is 


ft2  +  V(o*) 


q*  = 


k- 


(41) 


The  matrix  Q  is  a  diagonal  matrix  of  eigenfrequencies  for  the  perfectly  elastic  Earth;  q*  is  the 
kth  column  of  the  transformation  matrix  Q;  <7/,  is  the  complex  eigenfrequency  of  the  k,h  anelastic 
mode;  and  V ( <j)  is  an  anelastic  potential  energy  matrix  that  is  a  functional  of  products  of  the  elastic 
eigenfunctions  and  the  complex  frequency  dependent  elastic  moduli.  The  lateral  standing-wave 
nature  of  the  earth  free  oscillation  problem  leads  to  the  choice  of  the  frequency  as  the  dependent 
variable.  Because  the  elastic  moduli  depend  nonlinearly  on  frequency,  the  problem,  Eq.  (41),  for 
the  complex  eigenfrequencies  and  eigenvectors  is  nonlinear.  The  choice  of  wavenumber  as  the 
independent  variable  for  the  traveling  wave  problem  leads  to  the  linear  problem  we  have  derived 
above,  Eq.(40).  The  nonlinear  frequency  dependence  is  contained  in  the  elements  of  the  A,  B,  and 
C  matrices  of  Eq.(29),  Eq(30),  and  Eq.(31),  respectively. 

We  have  also  derived  similar  results  for  the  SH  problem.  We  form  the  composition 


dz{blZcn)  =  b[[MrE  +  SM]c, 


(42) 


As  definitions  of  b  and  M,  we  take  Eqs.  (16)  and  (30)  for  SH  motion.  For  M  we  take  Eq.  (17) 
with  k  replaced  by  K.  and  N  and  L  replaced  by  N  and  L,  respectively.  Likewise  the  definition  of  c 
follows  from  Eq.  (30)  and  (31).  Upon  integrating  Eq.  (42)  with  respect  to  z  from  0  to  —  °°,  carrying 
out  some  additional  algebra,  and  again  setting  co2  =  tJ2,  we  arrive  at 


(k2nl  -  k2N  -  L)  q„  =  K2  q„ 


(43) 
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where 


N 


5Nvi('Vm)rfe 


and 

dy  iW  dyi (m)  , 
dz  dz 


(44) 


(45) 


The  two  terms  <5L  and  5N  are  the  complex  frequency  dependent  parts  of  the  two  shear  moduli  L 
and  N.  The  anelasticity  tensor  c^/  for  a  linear  viscoelastic  material  can  be  written 


jki  —  jkl  +  8cijki((D).  (46) 

Substituting  the  appropriate  expressions  from  Eq.(46)  for  L  and  N,  we  were  able  to  separate  the 
perfectly  elastic  part  from  the  frequency  dependent  anelastic  part  for  the  SH  problem.  There  is 
no  restriction  on  the  magnitude  of  5L  and  5N,  making  it  possible  to  simplify  Eqs.(38)-(39).  It  is 
possible  to  also  write  the  complex  moduli  A,  C,  and  F  as  well  for  the  P-S  V  problem,  Eqs.(29)-(3 1), 
but  this  leads  to  algebraic  complexity  that  is  not  particularly  illuminating.  So  this  has  not  been 
done  for  the  P-SV  problem.  Also  note  that  Eq.  (43)  is  linear  in  K,,2.  A  final  point  is  that  both  Eq. 
(40)  and  Eq.  (43)  are  infinite  eigenvalue  equations.  Any  practical  implementation  will  necessarily 
employ  a  truncated  mode  set. 


4.  Discussion  and  Conclusions 

Advantages  of  using  the  elastic  eigenfunctions  as  a  basis  for  the  anelastic  eigenfunctions  are:  1. 
The  effect  of  anelasticity  on  individual  modes  can  be  examined  in  detail;  2.  Although  this  paper 
is  explicitly  for  laterally  homogenous  problems,  the  effect  of  range  dependent  attenuation  could 
be  studied  by  making  the  complex  expansion  coefficients  range  dependent.  (Pannatoni  (2011) 
has  published  a  coupled  mode  solution  for  a  range  dependent  all-fluid  acoustic  model,  which 
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includes  mode  coupling  due  to  attenuation.)  If  the  environment  is  not  geometrically  range  dependent 
(laterally  heterogeneous),  we  can  employ  the  same  elastic  basis;  3.  Used  in  conjunction  with  a 
general  range  dependent  coupled  mode  program,  the  propagation  physics  of  a  strongly  laterally 
heterogeneous  dependent  shallow  water  environment  or  regional  Earth  model  can  be  studied  in 
detail.  We  have  the  ability  to  isolate  the  influence  of  the  geometry,  and  different  aspects  of  the 
rheology  of  the  medium  on  a  propagating  seismic  or  seismo-acoustic  signal. 

We  derived  equations  Eq.  (40)  and  Eq.  (43)  for  the  computation  of  anelastic  surface-wave 
eigenfunctions  in  transversely  isotropic  media  by  expressing  them  in  terms  of  a  linear  complex 
generalized  eigenvalue  equation  for  the  P-SV  system,  and  a  linear  eigenvalue  equation  for  the  SH 
system.  No  perturbation  theory  is  needed.  The  completeness  of  the  elastic  modes  for  the  laterally 
homogeneous  Earth  as  a  basis  has  been  proved  by  Kirrmann  (1995)  for  the  case  of  a  free  surface  and 
rigid  lower  boundary,  which  is  the  locked  mode  approximation.  The  anelastic  modes  are  useful  for 
modeling  and  characterizing  seismo-acoustic  signals  propagating  in  a  shallow  water  environment 
or  regional  seismic  phases  characterized  by  high  attenuation  and  transverse  isotropy.  This  is  an 
environment  where  a  perturbation  treatment  of  the  bottom  or  upper  mantle  properties  applied  to 
mode  summation  signal  synthesis,  which  has  been  shown  by  previous  authors  to  lead  to  erroneous 
results.  The  solution  of  equations  (40)  and  (43)  are  used  to  represent  the  anelastic  modes  in  terms 
of  the  elastic  modes,  permitting  a  detailed  analysis  of  the  physics  of  strongly  bottom  interacting 
acoustic  propagation.  The  effects  of  transverse  isotropy  and  attenuation,  including  attenuation 
induced  dispersion,  are  properly  accounted  for.  Causality  is  assured  by  making  sure  the  complex, 
frequency  dependent  moduli  satisfy  the  Kramers-Kronig  relations.  Stable,  well-behaved  numerical 
algorithms  exist  for  solving  the  complex  generalized  eigenvalue  problem,  even  in  cases  where  the 
the  matrices  involved  are  near  singular  (Golub  and  Van  Loan,  1989).  The  QZ  method  is  suitable 
for  the  generalized  eigenvalue  problem.  Our  next  step  is  the  numerical  implementation  of  Eq.(40) 
and  Eq.(43).  A  suitable  check  would  be  a  camparison  with  the  results  obtained  from  Ivansson  and 
Karasalo(1992,  1993)  and  Ivansson’s  (1997)  direct  algorithm. 
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